A Trajectory Analysis using the ICESat2 Satellite

For information on adding the ICESat2 satellite to pygeodyn see the notebook titled, Adding a New Satellite: ICESat2.

Step 0: Pre-processing the PCE data into a G2B file

We will construct a G2B_PCEfile from a set of binary .rvg files which contain trajectory output from a reduced-dynamics run of GEODYN. The files will be stitched together to create a the GEODYN-specific, trajectory-based tracking data type (called PCE (Precice-Clock-Error) ) and stored in a G2B file.

Each .rvg file contains a 30-hour arc of ICESat2 trajectory data. These datasets are the output from a very precise run of GEODYN in which the orbit has been determined very well (to a few centimeters accuracy) using a reduced dynamics technique (empirical accelerations and other parameters are adjusted to account for any mismodelled forces).

The process that takes place in pygeodyn is as follows: 1. dump the data from each arc into a usable format 2. chop of the 3 hour padding on the ends to eliminate discontinuities from end effects 3. stitch together all the files 4. smooth over any existing discontinuities between arc gaps or maneuvers. 5. All data is placed in a single TRAJ.txt file which is then fed into a Fortran script (pce_converer.f) which converts the data into a G2B file to be ingested on fort.40 of the GEODYN run.

PygeodynPreprocessing

import copy

%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0,'/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/util_preprocessing/')
from PYGEODYN_Preprocess import PygeodynPreprocessing

path_to_prep_directory = '/data/data_geodyn/inputs/icesat2/pre_processing'
path_to_binaryrvg     = '/data/data_geodyn/inputs/icesat2/pre_processing/traj_files_rvg'

arc1_files = ['orbit.1807001.2018.312',
              'orbit.1807001.2018.313',
              'orbit.1807001.2018.314',
              'orbit.1807001.2018.315',
              'orbit.1807001.2018.316',
              'orbit.1807001.2018.317',
              'orbit.1807001.2018.318',
              'orbit.1807001.2018.319',
              'orbit.1807001.2018.320',
              'orbit.1807001.2018.321',
              'orbit.1807001.2018.322',
             ]

##### Uncomment the below call to overwrite the G2B data by re-running the pre-processing code
# Obj = PygeodynPreprocessing(path_to_binaryrvg, path_to_prep_directory,  arc1_files)
# Obj.run_preprocess_PCE()

Step 1: Run GEODYN using Pygeodyn with the ICESat2 configuration

As a reminder, PYGEODYN is called with the PYGEODYN class, but modifications to the ICESAT2 configuration are largely controlled through the Satellite_ICESat2 class in the PYGEODYN_Satellites.py file. Please refer to this file first for file structure setups and modifications to the Pygeodyn structure.

Step 1.1: Run with MSIS2

### Identify which arcs you want to run:
arcs_files = [
            '2018.313', # 1
            '2018.314', # 2
            '2018.315', # 3
            '2018.316', # 4
            '2018.317', # 5
#             '2018.318', # 6
#             '2018.319', # 7
#             '2018.320', # 8
#             '2018.321', # 9
#             '2018.322', # 10
             ]

#------ A dictionary containing the run parameters ------
run_params = {}
run_params['arc']              =   arcs_files
run_params['satellite']        =  'icesat2'
run_params['den_model']        =  'msis2'
run_params['SpecialRun_name']  =  '_TrajAnalysis'
run_params['verbose']          =  False
run_params['action']           =  'run'

%load_ext autoreload
%autoreload 2

sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/')
from PYGEODYN import Pygeodyn


#### ---------------------------------------------
####
#### ----------------- RUN MSIS2.0 ---------------
####
#### ---------------------------------------------

##### Use copy.deepcopy to copy all levels of dictionary and
###       allow modification of new variable
run_params1 = copy.deepcopy(run_params)
run_params1['den_model'] =  'msis2'

# ### Run pyeodyn for the arcs in the above set.
# Obj_Geodyn = Pygeodyn(run_params1)
# Obj_Geodyn.RUN_GEODYN()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Step 1.2: Run with DTM87 for comparison

%load_ext autoreload
%autoreload 2

sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/')
from PYGEODYN import Pygeodyn

#### ---------------------------------------------
####
#### ------------------- RUN DTM -----------------
####
#### ---------------------------------------------


##### Use copy.deepcopy to copy all levels of dictionary and
###       allow modification of new variable
run_params2 = copy.deepcopy(run_params)
run_params2['den_model']         =  'dtm87'

### Run pyeodyn for the arcs in the above set.
# Obj_Geodyn = Pygeodyn(run_params2)
# Obj_Geodyn.RUN_GEODYN()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Step 1.2: Run with Jaachia 71 Model for comparison

%load_ext autoreload
%autoreload 2

sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/')
from PYGEODYN import Pygeodyn


#### ----------------------------------------------
####
#### ----------------- RUN JAACHIA -----------------
####
#### ----------------------------------------------

##### Use copy.deepcopy to copy all levels of dictionary and
###       allow modification of new variable
run_params3 = copy.deepcopy(run_params)
run_params3['den_model']         =  'jaachia71'


##### Run pyeodyn for the arcs in the above set.
# Obj_Geodyn = Pygeodyn(run_params3)
# Obj_Geodyn.RUN_GEODYN()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Read GEODYN Output using PygeodynRead functionality

Get MSIS2 Data

# Clear the memory of the notebook because the saved variables get so large:
# %reset -f out

%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/')
from PYGEODYN import Pygeodyn

#------ A dictionary containing the run parameters ------
run_params1 = {}
run_params1['arc']               =  ['2018.313',
                                     '2018.314',
                                     '2018.315',
                                     '2018.316',
                                     '2018.317',
                                    ]
run_params1['satellite']         =  'icesat2'
run_params1['den_model']         =  'msis2'
run_params1['SpecialRun_name']   =  '_TrajAnalysis'
run_params1['verbose']           =  False
run_params1['action']            =  'read'

Obj_Geodyn1 = Pygeodyn(run_params1)
Obj_Geodyn1.getData()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
=================================================================
   Initializing PYGEODYN to ...
                            ... READ GEODYN output
=================================================================
     Loading ... icesat2_2018313_54hr.msis2
     Loading ... icesat2_2018314_54hr.msis2
     Loading ... icesat2_2018315_54hr.msis2
     Loading ... icesat2_2018316_54hr.msis2
     Loading ... icesat2_2018317_54hr.msis2

Get DTM87 Data

%load_ext autoreload
%autoreload 2

from PYGEODYN import Pygeodyn
import copy


run_params2 = copy.deepcopy(run_params1)

run_params2['den_model']         =  'dtm87'

Obj_Geodyn2 = Pygeodyn(run_params2)
Obj_Geodyn2.getData()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
=================================================================
   Initializing PYGEODYN to ...
                            ... READ GEODYN output
=================================================================
     Loading ... icesat2_2018313_54hr.dtm87
     Loading ... icesat2_2018314_54hr.dtm87
     Loading ... icesat2_2018315_54hr.dtm87
     Loading ... icesat2_2018316_54hr.dtm87
     Loading ... icesat2_2018317_54hr.dtm87

Get Jaachia 71 Data

%load_ext autoreload
%autoreload 2

from PYGEODYN import Pygeodyn
# import copy


run_params3 = copy.deepcopy(run_params1)
run_params3['den_model']         =  'jaachia71'

Obj_Geodyn3 = Pygeodyn(run_params3)
Obj_Geodyn3.getData()

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
=================================================================
   Initializing PYGEODYN to ...
                            ... READ GEODYN output
=================================================================
     Loading ... icesat2_2018313_54hr.jaachia71
     Loading ... icesat2_2018314_54hr.jaachia71
     Loading ... icesat2_2018315_54hr.jaachia71
     Loading ... icesat2_2018316_54hr.jaachia71
     Loading ... icesat2_2018317_54hr.jaachia71

Analyze output

Our analysis has the following output products:

1. Residuals of the POD across many arcs

2. RMS of fit of the POD across many arcs

3. Adjustment of the drag coefficient (drag acceleration to compensate for inaccuracies in the density model.)

4. Check of consistency (in the residuals) across overlapping arc times

5. RMS of the overlapping residual difference (this removes the PCE contribution) Using the Orbfil:

6. From ORBFIL grab the overlapping ephemeris and difference the two models. Compare this against the same of prediction

   • How well does the predicted time period match up with the determined ephemeris (see this in the resids of the two).

7. Calculate and plot the radial component of the trajectory

import plotly.graph_objects as go
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px


config = dict({
                'displayModeBar': True,
                'responsive': False,
#                 'staticPlot': True,
                'displaylogo': False,
                'showTips': False,
                })

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residual_meas_summary
from PYGEODYNAnalysis_icesat2PCEtrajectory import rms_summary_table


Obj_list = [Obj_Geodyn1,Obj_Geodyn2,Obj_Geodyn3,]
rms_summary_table(Obj_list)



fig = make_subplots(rows=2, cols=1,
     subplot_titles=(["Mean Residuals", 'RMS of Fit']),
     vertical_spacing = 0.1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn2, 0)
fig = plot_residual_meas_summary(fig, Obj_Geodyn3, 1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn1, 2)
fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
+———————————————————+—————————————————————————+————————————————+
|                     Summary Across all Arcs                  |
+———————————————————+—————————————————————————+————————————————+
+   Density Model   +   Mean Residual (cm)    +   RMS of Fit   +
+-------------------+-------------------------+----------------+
+       msis2       +       2.64000e+01       +  2.66840e-01   +
+       dtm87       +       1.14000e+01       +  5.02380e-01   +
+     jaachia71     +       2.20000e+00       +  2.61600e-01   +
+———————————————————+—————————————————————————+————————————————+

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residuals_observed



fig = make_subplots(rows=3, cols=1,
            subplot_titles=(['PCE X', 'PCE Y', 'PCE Z']),
            vertical_spacing = 0.1,
                       )

fig = plot_residuals_observed(fig, Obj_Geodyn2, 0)
fig = plot_residuals_observed(fig, Obj_Geodyn3, 1)
fig = plot_residuals_observed(fig, Obj_Geodyn1, 2)

fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

%load_ext autoreload
%autoreload 2
from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_cd_and_percdiff_from_apriori


fig = make_subplots(
    rows=2, cols=1,
    subplot_titles=(["Timeseries of Cd", "Percent difference from a priori (Cd=2.2)"]),
    vertical_spacing = 0.08,
    )
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn2, 0)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn3, 1)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn1, 2)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_ScaleDensity_with_CdScaleFactor__2

fig = make_subplots(rows=2, cols=1,
                subplot_titles=(["Model Ouptut Density", "Model Density * Cd Scaling Factor"]),
                shared_yaxes=True,
                vertical_spacing = 0.1,
                specs=[
                [{"secondary_y": False}],
                [{"secondary_y": False}], ])


fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn2, 0, 200)
fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn3, 1, 200)
fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn1, 2, 200)

# min_y = 1*1e-16
# max_y = 9*1e-12
# fig.update_yaxes(range=[min_y, max_y], row=1, col=1)
# fig.update_yaxes(range=[min_y, max_y], row=2, col=1)

fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Looking at the arc overlap time:

We want to show the residuals in the overlap time with the PCE data subtracted away.

####  ARC_OVERLAP_ObsResids_XYZ

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_XYZ

fig = make_subplots(rows=3, cols=1,
            subplot_titles=(['X', 'Y', 'Z']),
            vertical_spacing = 0.1,
            specs=[ [{"secondary_y": True }],
                    [{"secondary_y": True }],
                    [{"secondary_y": True }], ],)

arc1 = '2018.314'
arc2 = '2018.315'

fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn2, 0, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn3, 1, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn1, 2, arc1, arc2)

fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Residual Component Trajectory:

Convert the Interial XYZ coordinates to the satellite coordinate system (RSW). Then plot the radial component.

####  ARC_OVERLAP_ObsResids_XYZ

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_RSW_radial

fig = make_subplots(rows=2, cols=1,
            subplot_titles=(['Radial Component', 'Residual (PCE-ORBFIL)']),
            vertical_spacing = 0.2,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)

arc1 = '2018.314'
arc2 = '2018.315'

fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn2, 0, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn3, 1, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn1, 2, arc1, arc2)

fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload